{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "<a href=\"https://colab.research.google.com/github/CoreTheGreat/HBPU-Machine-Learning-Course/blob/main/ML_Chapter4_Clustering.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lPboLx_o0UxI"
   },
   "source": [
    "# 第四章：聚类\n",
    "湖北理工学院《机器学习》课程资料\n",
    "\n",
    "作者：李辉楚吴\n",
    "\n",
    "笔记内容概述:\n",
    "* 4.1 K均值聚类 k-Means\n",
    "* 4.2 密度聚类 DBSCAN\n",
    "* 4.3 电影评分分析中的聚类问题及应用\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Ifpm7Sql4U09"
   },
   "source": [
    "## 4.1 K均值聚类 k-Means\n",
    "\n",
    "本章需下载的数据：\n",
    "\n",
    "* 图片([下载](https://github.com/CoreTheGreat/HBPU-Machine-Learning-Course/tree/main/Data)): car_num.jpg, bmwk.png\n",
    "* MovieLens电影评分数据[下载](https://files.grouplens.org/datasets/movielens/ml-latest-small.zip)，具体说明见[此处](https://files.grouplens.org/datasets/movielens/ml-latest-small-README.html).\n",
    "\n",
    "截至2024年8越26日的官方介绍如下（结果可能与课堂讲解有出入）：\n",
    "These datasets will change over time, and are not appropriate for reporting research results. Small: 100,000 ratings and 3,600 tag applications applied to 9,000 movies by 600 users. Last updated 9/2018."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "sM-ziKb94S9_",
    "outputId": "a989e9b6-5ad2-4576-809f-90688d5fabd9"
   },
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "# Set the environment variable to avoid memory leak warning\n",
    "os.environ[\"OMP_NUM_THREADS\"] = \"2\"\n",
    "\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.cluster import DBSCAN\n",
    "from sklearn.datasets import make_moons\n",
    "from sklearn.datasets import make_blobs\n",
    "\n",
    "import time\n",
    "\n",
    "import cv2\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "color_list = ['tab:blue', 'tab:orange', 'tab:green', 'tab:red', 'tab:purple', 'tab:brown', 'tab:pink', 'tab:gray', 'tab:olive', 'tab:cyan']\n",
    "\n",
    "label_size = 18 # Label size\n",
    "ticklabel_size = 14 # Tick label size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UM_6PizbAnw3"
   },
   "source": [
    "### k-Means算法介绍：一维k-Means聚类"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 601
    },
    "id": "XuMonINmFXa7",
    "outputId": "25bbe732-d879-4b1b-9eeb-5139655f0474"
   },
   "outputs": [],
   "source": [
    "# Generate 1-D data\n",
    "np.random.seed(42)\n",
    "x = np.random.rand(100)\n",
    "\n",
    "# Plot\n",
    "fig, ax = plt.subplots(figsize=(7,7))\n",
    "ax.scatter(x, x, edgecolor='black', facecolor='white', linewidth=2, s=7**2)\n",
    "ax.tick_params(axis='both', which='major', labelsize=ticklabel_size) # Set tick label size\n",
    "\n",
    "# plt.savefig('1D_data_base.png', dpi=300) # Make figure clearer\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "386gEdEYGY4u",
    "outputId": "1f312c4d-7e38-47ea-8e13-4edb0e709cf7"
   },
   "outputs": [],
   "source": [
    "# k-Means clustering\n",
    "k = 2\n",
    "\n",
    "# Initialize cluster centers\n",
    "cluster_centers = [0.1, 0.2]\n",
    "\n",
    "# Initialize distance array\n",
    "distance = np.zeros((len(x), k))\n",
    "\n",
    "# Intercation\n",
    "max_iter = 10\n",
    "for iter_id in range(max_iter):\n",
    "    # Calculate distance between points to each center\n",
    "    for i in range(k):\n",
    "        distance[:,i] = np.abs(x - cluster_centers[i])\n",
    "\n",
    "    # Assign to closest centroid\n",
    "    cluster_idx = np.argmin(distance, axis=1)\n",
    "\n",
    "    # Display process\n",
    "    fig, ax = plt.subplots(figsize=(7,7))\n",
    "\n",
    "    # Display points close to center 1\n",
    "    ax.scatter(x[cluster_idx == 0], x[cluster_idx == 0], edgecolor=color_list[0], facecolor='white', linewidth=2, s=7**2)\n",
    "    ax.axvline(x=cluster_centers[0], color=color_list[0], linestyle='--', linewidth=2)\n",
    "\n",
    "    # Display points close to center 2\n",
    "    ax.scatter(x[cluster_idx == 1], x[cluster_idx == 1], edgecolor=color_list[1], facecolor='white', linewidth=2, s=7**2)\n",
    "    ax.axvline(x=cluster_centers[1], color=color_list[1], linestyle='--', linewidth=2)\n",
    "\n",
    "    ax.tick_params(axis='both', which='major', labelsize=ticklabel_size) # Set tick label size\n",
    "    ax.set_title(f'Interation: {iter_id}, C1: {cluster_centers[0]:.2f}, C2: {cluster_centers[1]:.2f}', fontsize=label_size)\n",
    "\n",
    "    # plt.savefig(f'1D_data_iter{iter_id}.png', dpi=300) # Make figure clearer\n",
    "    plt.show()\n",
    "\n",
    "    # Update cluster centers\n",
    "    for i in range(k):\n",
    "        cluster_centers[i] = np.mean(x[cluster_idx == i])\n",
    "\n",
    "# Display process\n",
    "fig, ax = plt.subplots(figsize=(7,7))\n",
    "\n",
    "# Display points close to center 1\n",
    "ax.scatter(x[cluster_idx == 0], x[cluster_idx == 0], edgecolor=color_list[0], facecolor='white', linewidth=2, s=7**2)\n",
    "ax.axvline(x=cluster_centers[0], color=color_list[0], linestyle='--', linewidth=2)\n",
    "\n",
    "# Display points close to center 2\n",
    "ax.scatter(x[cluster_idx == 1], x[cluster_idx == 1], edgecolor=color_list[1], facecolor='white', linewidth=2, s=7**2)\n",
    "ax.axvline(x=cluster_centers[1], color=color_list[1], linestyle='--', linewidth=2)\n",
    "\n",
    "ax.tick_params(axis='both', which='major', labelsize=ticklabel_size) # Set tick label size\n",
    "ax.set_title(f'Interation: {iter_id}, C1: {cluster_centers[0]:.2f}, C2: {cluster_centers[1]:.2f}', fontsize=label_size)\n",
    "\n",
    "# plt.savefig(f'1D_data_iter{iter_id}.png', dpi=300) # Make figure clearer\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "c-Xj4WtjPdwD"
   },
   "source": [
    "### 一维k-Means聚类的应用: 二值分割（Binary Segmentation）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 0
    },
    "id": "JWUUOgNeEyZ6",
    "outputId": "6d084a58-c228-4ed7-ec76-5d5a03f34679"
   },
   "outputs": [],
   "source": [
    "# Load the car number image\n",
    "image_carno = cv2.imread('./Data/car_num.jpg')\n",
    "\n",
    "# Convert from BGR to RGB\n",
    "image_carno = cv2.cvtColor(image_carno, cv2.COLOR_BGR2RGB)\n",
    "\n",
    "# Display the image\n",
    "fig, ax = plt.subplots(figsize=(7,7))\n",
    "img = ax.imshow(image_carno)\n",
    "plt.axis('off')  # Hide axes\n",
    "# plt.savefig('carno_base.png', dpi=300) # Make figure clearer\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 0
    },
    "id": "t8zeSOFyPs29",
    "outputId": "18c5d388-99d6-45a3-9715-9db7cf285116"
   },
   "outputs": [],
   "source": [
    "# Change image_carno into grey image\n",
    "image_carno_grey = cv2.cvtColor(image_carno, cv2.COLOR_BGR2GRAY)\n",
    "\n",
    "# Display the image\n",
    "fig, ax = plt.subplots(figsize=(7,7))\n",
    "img = ax.imshow(image_carno_grey, cmap='gray')\n",
    "plt.axis('off')  # Hide axes\n",
    "# plt.savefig('carno_grey.png', dpi=300) # Make figure clearer\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 0
    },
    "id": "MzSX7UoKQwyq",
    "outputId": "a5cf9244-4df4-4600-b435-1eda9b4adeba"
   },
   "outputs": [],
   "source": [
    "# Reshape image_carno_grey into a 1-D vector\n",
    "x = image_carno_grey.reshape(-1)\n",
    "\n",
    "# Using k-Means to separate background and foreground by pixels\n",
    "k = 2\n",
    "\n",
    "# Initialize cluster centers\n",
    "cluster_centers = np.random.rand(k) * 255\n",
    "start_centers = cluster_centers.copy()\n",
    "\n",
    "# Initialize distance array\n",
    "distance = np.zeros((len(x), k))\n",
    "\n",
    "# Intercation\n",
    "start_time = time.time()\n",
    "max_iter = 1000\n",
    "for iter_id in range(max_iter):\n",
    "    # Calculate distance between points to each center\n",
    "    for i in range(k):\n",
    "        distance[:,i] = np.abs(x - cluster_centers[i])\n",
    "\n",
    "    # Assign to closest centroid\n",
    "    cluster_idx = np.argmin(distance, axis=1)\n",
    "\n",
    "    # Update cluster centers\n",
    "    cluster_centers_prior = cluster_centers.copy()\n",
    "    for i in range(k):\n",
    "        cluster_centers[i] = np.mean(x[cluster_idx == i])\n",
    "\n",
    "    # Check if cluster_centers are stable enough to stop training\n",
    "    print(f'Iteration {iter_id}: Updated centers {cluster_centers}, Prior centers {cluster_centers_prior}')\n",
    "    if np.sum(np.abs(cluster_centers-cluster_centers_prior)) == 0:\n",
    "        break\n",
    "\n",
    "    cluster_centers_prior = cluster_centers\n",
    "\n",
    "end_time = time.time()\n",
    "print(f'Stop after iteration {iter_id}, time consumption is {end_time-start_time}')\n",
    "\n",
    "# Generate segmentation mask\n",
    "carno_mask_pixel = np.zeros_like(cluster_idx)\n",
    "low_value_cluster = np.argmin(cluster_centers)\n",
    "carno_mask_pixel[cluster_idx != low_value_cluster] = 1 # Set pixels with higher grey value to 1\n",
    "carno_mask_pixel = carno_mask_pixel.reshape(image_carno_grey.shape)\n",
    "\n",
    "# Display the mask\n",
    "fig, ax = plt.subplots(figsize=(7,7))\n",
    "img = ax.imshow(carno_mask_pixel, cmap='gray')\n",
    "ax.set_title(f'Final Centers: {cluster_centers}, Iteration: {iter_id}', fontsize=label_size)\n",
    "plt.axis('off')  # Hide axes\n",
    "# plt.savefig('carno_mask_pixel_2.png', dpi=300) # Make figure clearer\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "DJahqo4HcdkW"
   },
   "source": [
    "### 基于k-Means的图像压缩 (Image Compression)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 653
    },
    "id": "F0zcJGAle3fC",
    "outputId": "e27bc9f0-5afb-4039-fb64-8d57ab396f0c"
   },
   "outputs": [],
   "source": [
    "# Load the car number image\n",
    "image_bmwk = cv2.imread('./Data/bmwk.png')\n",
    "\n",
    "# Convert from BGR to RGB\n",
    "image_bmwk_rgb = cv2.cvtColor(image_bmwk, cv2.COLOR_BGR2RGB)\n",
    "\n",
    "x_r = image_bmwk_rgb[:, :, 0].reshape(-1) # Store colors in red channel\n",
    "x_g = image_bmwk_rgb[:, :, 1].reshape(-1) # Store colors in green channel\n",
    "x_b = image_bmwk_rgb[:, :, 2].reshape(-1) # Store colors in blue channel\n",
    "\n",
    "# Display the image with no margin\n",
    "plt.figure(figsize=(image_bmwk_rgb.shape[1]/100, image_bmwk_rgb.shape[0]/100))  # Convert pixels to inches\n",
    "plt.imshow(image_bmwk_rgb)\n",
    "plt.axis('off')  # Hide axes\n",
    "plt.subplots_adjust(left=0, right=1, top=1, bottom=0)  # Remove margins\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "Ya0qjoWejixW",
    "outputId": "77cf3d3d-02af-4af3-ba9d-8e5ae7d87354"
   },
   "outputs": [],
   "source": [
    "def KMeansImage(img, k):\n",
    "    # Get image size\n",
    "    w, h, c = img.shape\n",
    "\n",
    "    # Reshape image along channel\n",
    "    x = np.reshape(img, (w * h, c))\n",
    "\n",
    "    # Train k-Means model\n",
    "    mdl_km = KMeans(n_clusters=k, n_init='auto')\n",
    "    mdl_km.fit(x)\n",
    "\n",
    "    # Predict labels of each pixels\n",
    "    labels = mdl_km.predict(x).reshape(w, h)\n",
    "\n",
    "    # Get centers\n",
    "    center_colors = mdl_km.cluster_centers_ / 255.0\n",
    "\n",
    "    # Use center colors to generate compressed image\n",
    "    img_comp = np.zeros((w, h, c))\n",
    "    for i in range(w):\n",
    "        for j in range(h):\n",
    "            img_comp[i][j] = center_colors[labels[i][j]]\n",
    "    return img_comp, center_colors\n",
    "\n",
    "claster_num = [2, 4, 8, 16, 32, 64]\n",
    "for k in claster_num:\n",
    "    img_comp, center_colors = KMeansImage(image_bmwk_rgb, k)\n",
    "\n",
    "    # Display center colors\n",
    "    fig, ax = plt.subplots(figsize=(16,1))\n",
    "    ax.imshow([center_colors])\n",
    "    plt.axis('off')\n",
    "    # plt.savefig(f'bmwk_center_{k}.png', dpi=300)\n",
    "    plt.show()\n",
    "\n",
    "     # Display compressed image\n",
    "    plt.figure(figsize=(image_bmwk_rgb.shape[1]/100, image_bmwk_rgb.shape[0]/100))  # Convert pixels to inches\n",
    "    plt.imshow(img_comp)\n",
    "    plt.axis('off')\n",
    "    plt.subplots_adjust(left=0, right=1, top=1, bottom=0)  # Remove margins\n",
    "    # plt.savefig(f'bmwk_comp_{k}.png', format='png', bbox_inches='tight', pad_inches=0)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "DboOsOv-qBrg"
   },
   "source": [
    "### 手肘法确定k值"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 807
    },
    "id": "lwmavWf1xLkH",
    "outputId": "0d98d49f-750e-47b7-e868-35d7628bcaac"
   },
   "outputs": [],
   "source": [
    "# Using make_blobs to generate data of ten clustering\n",
    "X_mb, y_mb = make_blobs(n_samples=500, n_features=2, centers=6, random_state=42)\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "ax.scatter(X_mb[:, 0], X_mb[:, 1], marker=\"o\", c=y_mb, s=10**2, edgecolor=\"k\")\n",
    "plt.axis('off')\n",
    "# plt.savefig(f'make_blobs_base.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 552
    },
    "id": "ZDA08BSl08Bf",
    "outputId": "4e48a683-8890-41b5-87f7-26a87dabcf99"
   },
   "outputs": [],
   "source": [
    "k_list = np.arange(2, 20, 1)\n",
    "sse_list = np.zeros(len(k_list))\n",
    "\n",
    "mdl_km_list = []\n",
    "for i in range(len(k_list)):\n",
    "    mdl_km = KMeans(n_clusters=k_list[i], n_init='auto')\n",
    "    mdl_km.fit(X_mb)\n",
    "    mdl_km_list.append(mdl_km)\n",
    "    sse_list[i] = mdl_km.inertia_\n",
    "\n",
    "# Plot sse_list\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "ax.plot(k_list, sse_list, marker='o', linestyle='-', color='tab:blue')\n",
    "ax.set_xticks(k_list)\n",
    "ax.set_xlabel('Number of clusters (k)', fontsize=label_size)\n",
    "ax.set_ylabel('SSE', fontsize=label_size)\n",
    "ax.tick_params(axis='both', which='major', labelsize=ticklabel_size) # Set tick label size\n",
    "# plt.savefig(f'make_blobs_sse.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "K-JcsrtB23wb",
    "outputId": "1d4e1981-6201-4ef3-9752-907f558bd646"
   },
   "outputs": [],
   "source": [
    "# Display clustering result of k = 5, 6, 7\n",
    "k_disp = [5, 6, 7]\n",
    "for k in k_disp:\n",
    "    mdl_km = mdl_km_list[k-2]\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(10,10))\n",
    "    ax.scatter(X_mb[:, 0], X_mb[:, 1], marker=\"o\", c=mdl_km_list[k-2].labels_, s=10**2, edgecolor=\"k\")\n",
    "    ax.set_title(f'Number of clusters (k): {k}', fontsize=label_size)\n",
    "    plt.axis('off')\n",
    "\n",
    "    # plt.savefig(f'make_blobs_{k}.png', dpi=300)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Ky6e4lkilDdx"
   },
   "source": [
    "## 4.2 DBSCAN 聚类\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "k-Means聚类的问题，在Moon data上进行聚类"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 807
    },
    "id": "9qf2NtDdqFTl",
    "outputId": "d5552f56-4dc4-4d24-cd7b-0d81137e9955"
   },
   "outputs": [],
   "source": [
    "# Generate virtual data (Moon data)\n",
    "X_mm, y_mm = make_moons(n_samples=400, noise=0.05, random_state=42)\n",
    "\n",
    "# Normalize X_mm using z-score\n",
    "X_mm = (X_mm - np.min(X_mm, axis=0)) / (np.max(X_mm, axis=0)-np.min(X_mm, axis=0))\n",
    "\n",
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "ax.scatter(X_mm[:, 0], X_mm[:, 1], marker=\"o\", c=y_mm, s=10**2, edgecolor=\"k\")\n",
    "plt.axis('off')\n",
    "# plt.savefig(f'make_moon_base.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3xKX1XHnuXF2"
   },
   "source": [
    "用肘部法确定k值\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 552
    },
    "id": "4INr3adN6iN5",
    "outputId": "65351e31-e708-42e7-b500-78060475ad43"
   },
   "outputs": [],
   "source": [
    "k_list = np.arange(2, 20, 1)\n",
    "sse_list = np.zeros(len(k_list))\n",
    "\n",
    "mdl_km_list = []\n",
    "for i in range(len(k_list)):\n",
    "    mdl_km = KMeans(n_clusters=k_list[i], n_init='auto')\n",
    "    mdl_km.fit(X_mm)\n",
    "    mdl_km_list.append(mdl_km)\n",
    "    sse_list[i] = mdl_km.inertia_\n",
    "\n",
    "# Plot sse_list\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "ax.plot(k_list, sse_list, marker='o', linestyle='-', color='tab:blue')\n",
    "ax.set_xticks(k_list)\n",
    "ax.set_xlabel('Number of clusters (k)', fontsize=label_size)\n",
    "ax.set_ylabel('SSE', fontsize=label_size)\n",
    "ax.tick_params(axis='both', which='major', labelsize=ticklabel_size) # Set tick label size\n",
    "# plt.savefig(f'make_moon_sse.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "Y7wZOZqf698Y",
    "outputId": "ce51685e-cfa6-4aa8-af13-2ae76bc561c1"
   },
   "outputs": [],
   "source": [
    "# Display clustering result of k = 5, 6, 7\n",
    "k_disp = [2, 11, 19]\n",
    "for k in k_disp:\n",
    "    mdl_km = mdl_km_list[k-2]\n",
    "\n",
    "    fig, ax = plt.subplots(figsize=(10,10))\n",
    "    ax.scatter(X_mm[:, 0], X_mm[:, 1], marker=\"o\", c=mdl_km_list[k-2].labels_, s=10**2, edgecolor=\"k\")\n",
    "    ax.set_title(f'Number of clusters (k): {k}', fontsize=label_size)\n",
    "    plt.axis('off')\n",
    "\n",
    "    # plt.savefig(f'make_moon_{k}.png', dpi=300)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "noyzyu_yp16b"
   },
   "source": [
    "### DBSCAN 聚类基本逻辑"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-JIe2Q1XgTh2"
   },
   "source": [
    "数据分析与参数选择"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 551
    },
    "id": "g1hqSDWGgWah",
    "outputId": "af1c2d78-a4b6-4878-8d13-bb1bf4412f5e"
   },
   "outputs": [],
   "source": [
    "min_pts = 3\n",
    "\n",
    "# Comput Euclidean distance between samples\n",
    "distance_map = np.zeros((X_mm.shape[0], X_mm.shape[0]))\n",
    "for i in range(distance_map.shape[0]):\n",
    "    for j in range(distance_map.shape[1]):\n",
    "        distance_map[i,j] = np.linalg.norm(X_mm[i] - X_mm[j])\n",
    "\n",
    "# k-Distance\n",
    "k_distance = np.sort(np.sort(distance_map, axis=1)[:, min_pts])[::-1]\n",
    "\n",
    "# Draw k-Distance figure\n",
    "fig, ax = plt.subplots(figsize=(10,6))\n",
    "ax.plot(np.arange(1, len(k_distance)+1), k_distance, marker='o', linestyle='-', color='tab:blue')\n",
    "ax.set_xlabel('Sample Index', fontsize=label_size)\n",
    "ax.set_ylabel('k-Distance', fontsize=label_size)\n",
    "ax.tick_params(axis='both', which='major', labelsize=ticklabel_size) # Set tick label size\n",
    "# plt.savefig(f'dbscan_kdistance.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "gVhWKKp8xpRB"
   },
   "source": [
    "确定所有Core Point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "1t2vSctG8fQ2",
    "outputId": "d72dbb0f-6642-4399-9207-554c47ee0b69"
   },
   "outputs": [],
   "source": [
    "r_eps = 0.04\n",
    "\n",
    "cluster_labels = np.zeros(X_mm.shape[0]) - 1 # -1 means noise or unlabeled\n",
    "\n",
    "# Find all core points\n",
    "snap_noise_id = 0\n",
    "snap_corepoint_id = 0\n",
    "for i in range(len(cluster_labels)):\n",
    "    if np.sum(distance_map[i, :] <= r_eps) >= min_pts:\n",
    "        cluster_labels[i] = 0 # 0 means unclustered core point\n",
    "\n",
    "        if snap_corepoint_id < 3:\n",
    "            fig, ax = plt.subplots(figsize=(10,10))\n",
    "\n",
    "            # Draw unlabeled or noise point\n",
    "            noise_idx = (cluster_labels == -1)\n",
    "            ax.scatter(X_mm[noise_idx, 0], X_mm[noise_idx, 1], marker=\"o\", c=\"k\", s=10**2, edgecolor=\"k\", zorder=0)\n",
    "\n",
    "            # Draw eps field of the current core point\n",
    "            corepoint_idx = (cluster_labels == 0)\n",
    "            circle = plt.Circle((X_mm[i, 0], X_mm[i, 1]), r_eps, edgecolor='red', facecolor='tab:red', alpha=0.5, zorder=1)\n",
    "            ax.add_patch(circle)\n",
    "\n",
    "            # Draw core points\n",
    "            ax.scatter(X_mm[corepoint_idx, 0], X_mm[corepoint_idx, 1], marker=\"o\", c='red', s=10**2, edgecolor=\"k\", zorder=2)\n",
    "\n",
    "            ax.set_title(f'Phase 1: Iteration {i}, eps = {r_eps}, minPts = {min_pts}', fontsize=label_size)\n",
    "            plt.axis('off')\n",
    "            # plt.savefig(f'dbscan_corepoint_{snap_corepoint_id+1}.png', dpi=300)\n",
    "            plt.show()\n",
    "\n",
    "            snap_corepoint_id += 1\n",
    "    elif snap_noise_id < 3:\n",
    "        fig, ax = plt.subplots(figsize=(10,10))\n",
    "\n",
    "        # Draw unlabeled or noise point\n",
    "        noise_idx = (cluster_labels == -1)\n",
    "        ax.scatter(X_mm[noise_idx, 0], X_mm[noise_idx, 1], marker=\"o\", c=\"k\", s=10**2, edgecolor=\"k\", zorder=0)\n",
    "\n",
    "        # Draw core points\n",
    "        corepoint_idx = (cluster_labels == 0)\n",
    "        ax.scatter(X_mm[corepoint_idx, 0], X_mm[corepoint_idx, 1], marker=\"o\", c='red', s=10**2, edgecolor=\"k\", zorder=2)\n",
    "\n",
    "        # Draw eps field of the current point\n",
    "        circle = plt.Circle((X_mm[i, 0], X_mm[i, 1]), r_eps, edgecolor='blue', facecolor='tab:blue', alpha=0.5, zorder=1)\n",
    "        ax.add_patch(circle)\n",
    "\n",
    "        ax.set_title(f'Phase 1: Iteration {i}, eps = {r_eps}, minPts = {min_pts}', fontsize=label_size)\n",
    "        plt.axis('off')\n",
    "        # plt.savefig(f'dbscan_noise_{snap_noise_id+1}.png', dpi=300)\n",
    "        plt.show()\n",
    "\n",
    "        snap_noise_id += 1\n",
    "\n",
    "# Cluster all core points\n",
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "\n",
    "# Draw unlabeled or noise point\n",
    "noise_idx = (cluster_labels == -1)\n",
    "ax.scatter(X_mm[noise_idx, 0], X_mm[noise_idx, 1], marker=\"o\", c=\"k\", s=10**2, edgecolor=\"k\", zorder=0)\n",
    "\n",
    "# Draw core points\n",
    "corepoint_idx = (cluster_labels == 0)\n",
    "ax.scatter(X_mm[corepoint_idx, 0], X_mm[corepoint_idx, 1], marker=\"o\", c='red', s=10**2, edgecolor=\"k\", zorder=2)\n",
    "\n",
    "ax.set_title(f'Phase 1: Final, eps = {r_eps}, minPts = {min_pts}', fontsize=label_size)\n",
    "plt.axis('off')\n",
    "# plt.savefig(f'dbscan_corepoint_final.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "y05jO6mnxw4h"
   },
   "source": [
    "将互相可达的核心点合并成簇"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "olhR5Ow4yLEP",
    "outputId": "4b04e102-7722-4544-aabb-261be4ad7725"
   },
   "outputs": [],
   "source": [
    "# Scan unlabeled core points\n",
    "corepoint_cluster_labels = cluster_labels.copy()\n",
    "corepoint_indices = np.where(corepoint_cluster_labels == 0)[0]\n",
    "\n",
    "snap_corepoint_cluster_itr_idx = 0\n",
    "snap_corepoint_cluster_idx = 0\n",
    "\n",
    "cluster_id = 0 # Init cluster label\n",
    "\n",
    "# Clustering until all core points were successfully asigned labels\n",
    "while corepoint_cluster_labels[corepoint_indices].min() == 0:\n",
    "\n",
    "    cluster_id += 1 # Build new cluster\n",
    "\n",
    "    # Find start point, each cluster will be built completed in each loop\n",
    "    for start_corepoint_idx in corepoint_indices:\n",
    "        if corepoint_cluster_labels[start_corepoint_idx] == 0:\n",
    "            candidate_idx = np.array([start_corepoint_idx]) # Regard the single point as a cluster\n",
    "            break\n",
    "\n",
    "    # Repeat until no unlabeled points were found\n",
    "    while candidate_idx.size > 0:\n",
    "        if cluster_id == 2 and snap_corepoint_cluster_idx < 1:\n",
    "            fig, ax = plt.subplots(figsize=(10,10))\n",
    "\n",
    "            # Draw noise points\n",
    "            noise_idx = (corepoint_cluster_labels == -1)\n",
    "            ax.scatter(X_mm[noise_idx, 0], X_mm[noise_idx, 1], marker=\"o\", c=\"k\", s=10**2, edgecolor=\"k\", zorder=0)\n",
    "\n",
    "            # Draw unlabeled core points\n",
    "            unlabel_idx = (corepoint_cluster_labels == 0)\n",
    "            ax.scatter(X_mm[unlabel_idx, 0], X_mm[unlabel_idx, 1], marker=\"o\", c=\"r\", s=10**2, edgecolor=\"k\", zorder=0)\n",
    "\n",
    "            # Draw cluster 1 core points\n",
    "            cluster1_idx = (corepoint_cluster_labels == 1)\n",
    "            ax.scatter(X_mm[cluster1_idx, 0], X_mm[cluster1_idx, 1], marker=\"o\", c='tab:blue', s=10**2, edgecolor=\"k\", zorder=2)\n",
    "\n",
    "            # Draw cluster 2 core points\n",
    "            cluster2_idx = (corepoint_cluster_labels == 2)\n",
    "            ax.scatter(X_mm[cluster2_idx, 0], X_mm[cluster2_idx, 1], marker=\"o\", c='yellow', s=10**2, edgecolor=\"k\", zorder=2)\n",
    "\n",
    "            ax.set_title(f'Phase 2: Cluster 1, eps = {r_eps}, minPts = {min_pts}', fontsize=label_size)\n",
    "            plt.axis('off')\n",
    "            # plt.savefig(f'dbscan_corepoint_cluster_1.png', dpi=300)\n",
    "            plt.show()\n",
    "\n",
    "            snap_corepoint_cluster_idx += 1\n",
    "\n",
    "        # Assign cluster label to the first candidate core point\n",
    "        corepoint_cluster_labels[candidate_idx[0]] = cluster_id\n",
    "        circle_id = candidate_idx[0]\n",
    "\n",
    "        # Find the neighbor core points of the first candidate core point\n",
    "        neighbor_idx = np.where(distance_map[candidate_idx[0]] <= r_eps)[0] # Find neighbor points\n",
    "        neighbor_idx = neighbor_idx[neighbor_idx != candidate_idx[0]] # Filter out itself\n",
    "        neighbor_idx = neighbor_idx[corepoint_cluster_labels[neighbor_idx] == 0] # Remain unlabeled core points\n",
    "\n",
    "        # Concatenate neighbor_idx to candidate_idx\n",
    "        candidate_idx = np.concatenate((candidate_idx, neighbor_idx))\n",
    "\n",
    "        # Remove the first candidate_idx core point\n",
    "        candidate_idx = np.delete(candidate_idx, 0)\n",
    "\n",
    "        # Remove duplicated points\n",
    "        candidate_idx = np.unique(candidate_idx)\n",
    "\n",
    "        if snap_corepoint_cluster_itr_idx < 3:\n",
    "            fig, ax = plt.subplots(figsize=(10,10))\n",
    "\n",
    "            # Draw noise points\n",
    "            noise_idx = (corepoint_cluster_labels == -1)\n",
    "            ax.scatter(X_mm[noise_idx, 0], X_mm[noise_idx, 1], marker=\"o\", c=\"k\", s=10**2, edgecolor=\"k\", zorder=0)\n",
    "\n",
    "            # Draw unlabeled core points\n",
    "            unlabel_idx = (corepoint_cluster_labels == 0)\n",
    "            ax.scatter(X_mm[unlabel_idx, 0], X_mm[unlabel_idx, 1], marker=\"o\", c=\"r\", s=10**2, edgecolor=\"k\", zorder=0)\n",
    "\n",
    "            # Draw eps field of the current core point\n",
    "            circle = plt.Circle((X_mm[circle_id, 0], X_mm[circle_id, 1]), r_eps, edgecolor='blue', facecolor='tab:blue', alpha=0.5, zorder=1)\n",
    "            ax.add_patch(circle)\n",
    "\n",
    "            # Draw cluster 1 core points\n",
    "            cluster1_idx = (corepoint_cluster_labels == 1)\n",
    "            ax.scatter(X_mm[cluster1_idx, 0], X_mm[cluster1_idx, 1], marker=\"o\", c='tab:blue', s=10**2, edgecolor=\"k\", zorder=2)\n",
    "\n",
    "            # Draw cluster 2 core points\n",
    "            cluster2_idx = (corepoint_cluster_labels == 2)\n",
    "            ax.scatter(X_mm[cluster2_idx, 0], X_mm[cluster2_idx, 1], marker=\"o\", c='yellow', s=10**2, edgecolor=\"k\", zorder=2)\n",
    "\n",
    "            ax.set_title(f'Phase 2: Iteration {snap_corepoint_cluster_itr_idx+1}, eps = {r_eps}, minPts = {min_pts}', fontsize=label_size)\n",
    "            plt.axis('off')\n",
    "            # plt.savefig(f'dbscan_corepoint_cluster_Itr{snap_corepoint_cluster_itr_idx}.png', dpi=300)\n",
    "            plt.show()\n",
    "\n",
    "            snap_corepoint_cluster_itr_idx += 1\n",
    "\n",
    "# Display clustering result\n",
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "\n",
    "# Draw noise points\n",
    "noise_idx = (corepoint_cluster_labels == -1)\n",
    "ax.scatter(X_mm[noise_idx, 0], X_mm[noise_idx, 1], marker=\"o\", c=\"k\", s=10**2, edgecolor=\"k\", zorder=0)\n",
    "\n",
    "# Draw unlabeled core points\n",
    "unlabel_idx = (corepoint_cluster_labels == 0)\n",
    "ax.scatter(X_mm[unlabel_idx, 0], X_mm[unlabel_idx, 1], marker=\"o\", c=\"r\", s=10**2, edgecolor=\"k\", zorder=0)\n",
    "\n",
    "# Draw cluster 1 core points\n",
    "cluster1_idx = (corepoint_cluster_labels == 1)\n",
    "ax.scatter(X_mm[cluster1_idx, 0], X_mm[cluster1_idx, 1], marker=\"o\", c='tab:blue', s=10**2, edgecolor=\"k\", zorder=2)\n",
    "\n",
    "# Draw cluster 2 core points\n",
    "cluster2_idx = (corepoint_cluster_labels >= 2)\n",
    "ax.scatter(X_mm[cluster2_idx, 0], X_mm[cluster2_idx, 1], marker=\"o\", c=corepoint_cluster_labels[cluster2_idx], s=10**2, edgecolor=\"k\", zorder=2)\n",
    "\n",
    "ax.set_title(f'Phase 2: Final, eps = {r_eps}, minPts = {min_pts}', fontsize=label_size)\n",
    "plt.axis('off')\n",
    "# plt.savefig(f'dbscan_corepoint_cluster.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QvcXFixnVxyJ"
   },
   "source": [
    "将满足条件的边缘点Border Point加入簇中"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "3KS5UbePo-_Z",
    "outputId": "7c573939-e110-4e52-9325-f6003f255f75"
   },
   "outputs": [],
   "source": [
    "borderpoint_labels = corepoint_cluster_labels.copy()\n",
    "\n",
    "# Get unlabeled points\n",
    "candidate_indices = np.where(borderpoint_labels == -1)[0]\n",
    "\n",
    "snap_border_itr = 0\n",
    "\n",
    "# Scan candidate_idx to find border points\n",
    "for candidate_idx in candidate_indices:\n",
    "\n",
    "    # Find the neighbor points of i\n",
    "    neighbor_idx = np.where(distance_map[candidate_idx] <= r_eps)[0] # Find neighbor points\n",
    "    neighbor_idx = neighbor_idx[neighbor_idx != candidate_idx] # Filter out itself\n",
    "    neighbor_idx = neighbor_idx[borderpoint_labels[neighbor_idx] > 0] # Remain core points\n",
    "\n",
    "    # Candidate point is a border point\n",
    "    if len(neighbor_idx) > 0:\n",
    "\n",
    "        # Assign the closest core point's cluster label to the candicate point\n",
    "        distance_to_corepoint = distance_map[candidate_idx, neighbor_idx]\n",
    "        closest_corepoint_idx = neighbor_idx[np.argmin(distance_to_corepoint)]\n",
    "        borderpoint_labels[candidate_idx] = borderpoint_labels[closest_corepoint_idx]\n",
    "\n",
    "        if snap_border_itr < 3:\n",
    "            # Display clustering result\n",
    "            fig, ax = plt.subplots(figsize=(10,10))\n",
    "\n",
    "            # Draw noise points\n",
    "            noise_idx = (borderpoint_labels == -1)\n",
    "            ax.scatter(X_mm[noise_idx, 0], X_mm[noise_idx, 1], marker=\"o\", c=\"k\", s=10**2, edgecolor=\"k\", zorder=0)\n",
    "\n",
    "            # Draw eps field of the current core point\n",
    "            circle = plt.Circle((X_mm[candidate_idx, 0], X_mm[candidate_idx, 1]), r_eps, edgecolor='red', facecolor='tab:red', alpha=0.5, zorder=1)\n",
    "            ax.add_patch(circle)\n",
    "\n",
    "            # Draw cluster 1 core points\n",
    "            cluster1_idx = (borderpoint_labels == 1)\n",
    "            ax.scatter(X_mm[cluster1_idx, 0], X_mm[cluster1_idx, 1], marker=\"o\", c='tab:blue', s=10**2, edgecolor=\"k\", zorder=2)\n",
    "\n",
    "            # Draw cluster 2 core points\n",
    "            cluster2_idx = (borderpoint_labels >= 2)\n",
    "            ax.scatter(X_mm[cluster2_idx, 0], X_mm[cluster2_idx, 1], marker=\"o\", c=borderpoint_labels[cluster2_idx], s=10**2, edgecolor=\"k\", zorder=2)\n",
    "\n",
    "            ax.set_title(f'Phase 3: Itr {snap_border_itr}, eps = {r_eps}, minPts = {min_pts}', fontsize=label_size)\n",
    "            plt.axis('off')\n",
    "            # plt.savefig(f'dbscan_border_{snap_border_itr}.png', dpi=300)\n",
    "            plt.show()\n",
    "\n",
    "            snap_border_itr += 1\n",
    "\n",
    "# Display clustering result\n",
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "\n",
    "# Draw noise points\n",
    "noise_idx = (borderpoint_labels == -1)\n",
    "ax.scatter(X_mm[noise_idx, 0], X_mm[noise_idx, 1], marker=\"o\", c=\"k\", s=10**2, edgecolor=\"k\", zorder=0)\n",
    "\n",
    "# Draw cluster 1 core points\n",
    "cluster1_idx = (borderpoint_labels == 1)\n",
    "ax.scatter(X_mm[cluster1_idx, 0], X_mm[cluster1_idx, 1], marker=\"o\", c='tab:blue', s=10**2, edgecolor=\"k\", zorder=2)\n",
    "\n",
    "# Draw cluster 2 core points\n",
    "cluster2_idx = (borderpoint_labels >= 2)\n",
    "ax.scatter(X_mm[cluster2_idx, 0], X_mm[cluster2_idx, 1], marker=\"o\", c=borderpoint_labels[cluster2_idx], s=10**2, edgecolor=\"k\", zorder=2)\n",
    "\n",
    "ax.set_title(f'Phase 3: Final, eps = {r_eps}, minPts = {min_pts}', fontsize=label_size)\n",
    "plt.axis('off')\n",
    "# plt.savefig(f'dbscan_border_final.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EPJRDqroj_na"
   },
   "source": [
    "人工调整参数后的结果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 850
    },
    "id": "9qRRN2OBkDl5",
    "outputId": "5f6e35b8-362c-46fd-d537-4babfee454e1"
   },
   "outputs": [],
   "source": [
    "mdl_dbscan_eps = 0.05\n",
    "mdl_dbscan_minpts = 3\n",
    "\n",
    "# Define DBSCAN model\n",
    "mdl_dbscan = DBSCAN(eps=mdl_dbscan_eps, min_samples=mdl_dbscan_minpts)\n",
    "\n",
    "# Train model\n",
    "start_time = time.time()\n",
    "mdl_dbscan.fit(X_mm)\n",
    "end_time = time.time()\n",
    "\n",
    "print(f'Training time: {end_time - start_time:.2f} seconds')\n",
    "\n",
    "labels_mdl_dbscan = mdl_dbscan.labels_\n",
    "\n",
    "# Display clustering result\n",
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "\n",
    "# Draw noise points\n",
    "noise_idx = (labels_mdl_dbscan == -1)\n",
    "ax.scatter(X_mm[noise_idx, 0], X_mm[noise_idx, 1], marker=\"o\", c=\"k\", s=10**2, edgecolor=\"k\", zorder=0)\n",
    "\n",
    "# Draw cluster 1 core points\n",
    "cluster1_idx = (labels_mdl_dbscan == 0)\n",
    "ax.scatter(X_mm[cluster1_idx, 0], X_mm[cluster1_idx, 1], marker=\"o\", c='tab:blue', s=10**2, edgecolor=\"k\", zorder=2)\n",
    "\n",
    "# Draw cluster 2 core points\n",
    "cluster2_idx = (labels_mdl_dbscan >= 1)\n",
    "ax.scatter(X_mm[cluster2_idx, 0], X_mm[cluster2_idx, 1], marker=\"o\", c=labels_mdl_dbscan[cluster2_idx], s=10**2, edgecolor=\"k\", zorder=2)\n",
    "\n",
    "ax.set_title(f'Adjust parameter manually: eps = {mdl_dbscan_eps}, minPts = {mdl_dbscan_minpts}', fontsize=label_size)\n",
    "plt.axis('off')\n",
    "# plt.savefig(f'mdl_dbscan.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Kt6Fk7eTyAXW"
   },
   "source": [
    "## 4.3 电影评分分析中的聚类问题及应用"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "WYv0kOOZWEaT"
   },
   "source": [
    "Moves Data File Structure (movies.csv)\n",
    "\n",
    "Genres: (no genres listed), 动作（Action）, 冒险（Adventure）, 动画（Animation）, 儿童（Children）, 喜剧（Comedy）, 犯罪（Crime）, 纪录片（Documentary）, 剧情（Drama）, 奇幻（Fantasy）, 黑色电影（Film-Noir）, 恐怖（Horror）, IMAX, 音乐（Musical）, 推理（Mystery）, 爱情（Romance）, 科幻（Sci-Fi）, 惊悚（Thriller）, 战争（War）, 西部（Western）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "1JbjUpmKUhuv",
    "outputId": "0fb8b51a-97ed-4855-cb89-71528b2d08ed"
   },
   "outputs": [],
   "source": [
    "movies = pd.read_csv('./Data/movies.csv')\n",
    "print(f'Movie number {len(movies)}')\n",
    "movies.head(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Statistic genres\n",
    "genres_list = []\n",
    "for genres in movies['genres']:\n",
    "    genres_list.extend(genres.split('|'))\n",
    "genres_list = np.sort(np.unique(genres_list))\n",
    "print(f'Movie genres: {genres_list}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "XLSroDN4Y7CY"
   },
   "source": [
    "Rating Data File Structure (ratings.csv)\n",
    "\n",
    "Each line of this file after the header row represents one rating of one movie by one user, and has the following format:\n",
    "\n",
    "userId, movieId, rating, timestamp\n",
    "\n",
    "Ratings are made on a 5-star scale, with half-star increments (0.5 stars - 5.0 stars).\n",
    "\n",
    "Timestamps represent seconds since midnight Coordinated Universal Time (UTC) of January 1, 1970."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 224
    },
    "id": "05vhvNwyUujv",
    "outputId": "3c3085a2-7acc-4f0c-df1c-7b343ff71aae"
   },
   "outputs": [],
   "source": [
    "ratings = pd.read_csv('./Data/ratings.csv')\n",
    "user_list = np.sort(np.unique(ratings['userId']))\n",
    "print(f'{len(user_list)} have provided {len(ratings)} rate records')\n",
    "ratings.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "fdYsKkr4gD1d"
   },
   "source": [
    "统计用户对各种电影类型的喜爱程度，即计算用户对各类电影的平均分"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "kK7lw-OTEqM3",
    "outputId": "508f6d21-f647-4ab0-a497-12cd075d7330"
   },
   "outputs": [],
   "source": [
    "# Creat user-genres array\n",
    "user_genres_rate = np.zeros((len(user_list), len(genres_list)))\n",
    "user_genres_rate_counts = np.zeros((len(user_list), len(genres_list)))\n",
    "\n",
    "for i in range(len(ratings)):\n",
    "    # User ID start from 1\n",
    "    # Convert user ID to user_idx by minus 1\n",
    "    user_idx = int(ratings.iloc[i]['userId'] - 1)\n",
    "    \n",
    "    # Get rate\n",
    "    rate = ratings.iloc[i]['rating']\n",
    "\n",
    "    # Get target movie\n",
    "    movie_id = ratings.iloc[i]['movieId']\n",
    "    \n",
    "    # Split movie genres string\n",
    "    genres_type = movies[movies['movieId'] == movie_id]['genres'].values[0]\n",
    "    genres_type = genres_type.split('|')\n",
    "\n",
    "    # Statistic rates of movie genres\n",
    "    for genres in genres_type:\n",
    "        # Using '[0][0]' to get index from a list tuple\n",
    "        # First [0] get the index list outputed by np.where()\n",
    "        # Second [0] get the first item of list (with single item)\n",
    "        genres_idx = np.where(genres_list == genres)[0][0]\n",
    "\n",
    "        # Sum rate\n",
    "        user_genres_rate[user_idx, genres_idx] += rate\n",
    "        \n",
    "        # Count movie number of the genres\n",
    "        user_genres_rate_counts[user_idx, genres_idx] += 1\n",
    "\n",
    "# Compute genres rates of all users\n",
    "user_genres_rate = user_genres_rate / user_genres_rate_counts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "观察user_genres_rate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Display user_genres_rate shape in a table\n",
    "import pandas as pd\n",
    "\n",
    "# Convert user_genres_rate to a DataFrame\n",
    "user_genres_df = pd.DataFrame(user_genres_rate, columns=genres_list)\n",
    "\n",
    "# Display shape information\n",
    "print(\"\\nDataFrame Shape:\")\n",
    "print(f\"Rows: {user_genres_df.shape[0]}\")\n",
    "print(f\"Columns: {user_genres_df.shape[1]}\")\n",
    "\n",
    "# Display the first few rows of the DataFrame, showing only 'Action' and 'Animation'\n",
    "user_genres_df[['Action', 'Animation']].head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "用散点图展示动作电影（Action）和动画电影（Animation）的用户喜爱情况\n",
    "\n",
    "可选电影类型: \n",
    "\n",
    "'(no genres listed)' 'Action' 'Adventure' 'Animation' 'Children' 'Comedy'\n",
    "\n",
    "'Crime' 'Documentary' 'Drama' 'Fantasy' 'Film-Noir' 'Horror' 'IMAX'\n",
    "\n",
    "'Musical' 'Mystery' 'Romance' 'Sci-Fi' 'Thriller' 'War' 'Western'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "KaY4gClVJB6v",
    "outputId": "8527cef2-c3bb-4ca0-e584-9e328a88d557"
   },
   "outputs": [],
   "source": [
    "genres_1 = 'Action' # Define first genres\n",
    "genres_2 = 'Animation' # Define second genres\n",
    "\n",
    "idx_g1 = np.where(genres_list == genres_1)[0][0] # Get genres 1 index in the genres_list\n",
    "idx_g2 = np.where(genres_list == genres_2)[0][0] # Get genres 2 index in the genres_list\n",
    "\n",
    "# Get list of user who have watched movies including both genres_1 and genres_2\n",
    "# Idea: People who have watched related movies are more credible in determining whether they like them or not.\n",
    "idx_both = np.where(np.logical_and(user_genres_rate[:,idx_g1] > 0, user_genres_rate[:,idx_g2] > 0))[0]\n",
    "print(f'Number of users who watched both {genres_1} and {genres_2}: {len(idx_both)}')\n",
    "\n",
    "# Filter rating information of users who have watched both genres_1 and genres_2\n",
    "x = user_genres_rate[idx_both][:, [idx_g1, idx_g2]]\n",
    "print(f'Shape of x: {x.shape}')\n",
    "\n",
    "# Drawing rates distribution by scatter figure\n",
    "fig, ax = plt.subplots(figsize=(6,6))\n",
    "\n",
    "ax.scatter(x[:,0], x[:,1], marker=\"o\", c='white', s=10**2, edgecolor=\"k\")\n",
    "\n",
    "ax.tick_params(axis='both', which='major', labelsize=ticklabel_size) # Set tick label size\n",
    "ax.set_xlabel(f'Average rate of {genres_1}', fontsize=label_size)\n",
    "ax.set_ylabel(f'Average rate of {genres_2}', fontsize=label_size)\n",
    "ax.set_xlim([-0.1, 5.1])\n",
    "ax.set_ylim([-0.1, 5.1])\n",
    "ax.tick_params(axis='both', which='major', labelsize=ticklabel_size) # Set tick label size\n",
    "# plt.savefig(f'rate_of_interested_movies.png', dpi=300) # Make figure clearer\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用K-Means分类，观察结果，进行分析"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for k in [2, 4, 8]:\n",
    "\n",
    "    # Define k-means clustering model\n",
    "    mdl_km = KMeans(n_clusters=k, n_init='auto')\n",
    "\n",
    "    # Train k-Means model\n",
    "    mdl_km.fit(x)\n",
    "    \n",
    "    # Assign clusters\n",
    "    labels_km = mdl_km.predict(x)\n",
    "\n",
    "    # Draw scatter figure\n",
    "    fig, ax = plt.subplots(figsize=(6,6))\n",
    "\n",
    "    ax.scatter(x[:,0], x[:,1], marker=\"o\", c=labels_km, s=10**2, edgecolor=\"k\")\n",
    "\n",
    "    ax.tick_params(axis='both', which='major', labelsize=ticklabel_size) # Set tick label size\n",
    "    ax.set_xlabel(f'Average rate of {genres_1}', fontsize=label_size)\n",
    "    ax.set_ylabel(f'Average rate of {genres_2}', fontsize=label_size)\n",
    "    ax.set_xlim([-0.1, 5.1])\n",
    "    ax.set_ylim([-0.1, 5.1])\n",
    "    ax.set_title(f'k-Means clustering, k = {k}', fontsize=label_size)\n",
    "    ax.tick_params(axis='both', which='major', labelsize=ticklabel_size) # Set tick label size\n",
    "    plt.savefig(f'rate_kMeans_{k}.png', dpi=300) # Make figure clearer\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用DBSCAN分类，观察结果，进行分析"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 551
    },
    "id": "TGSPq0ZlwU_V",
    "outputId": "851d99ed-fc4e-4768-e9ad-78e6c0b624a3"
   },
   "outputs": [],
   "source": [
    "# Using k-Distance to find potential values of e\n",
    "min_pts = 3\n",
    "\n",
    "# Initialize distance_map to store distance\n",
    "distance_map = np.zeros((len(x), len(x)))\n",
    "\n",
    "# Compute Manhattan distance using nested loops\n",
    "for i in range(len(x)):\n",
    "    for j in range(len(x)):\n",
    "        # Calculate Manhattan distance between point i and point j\n",
    "        distance_map[i, j] = np.abs(x[i, 0] - x[j, 0]) + np.abs(x[i, 1] - x[j, 1])\n",
    "\n",
    "# Get k-Distance from distance_map\n",
    "k_distance = np.sort(distance_map, axis=1)[:, min_pts]\n",
    "\n",
    "# Sort k_distance in descend order\n",
    "k_distance = np.sort(k_distance)[::-1]\n",
    "\n",
    "# Draw k-Distance figure\n",
    "fig, ax = plt.subplots(figsize=(8,6))\n",
    "ax.plot(np.arange(1, len(k_distance)+1), k_distance, marker='o', linestyle='-', color='tab:blue')\n",
    "ax.set_xlabel('Sample Index', fontsize=label_size)\n",
    "ax.set_ylabel('k-Distance', fontsize=label_size)\n",
    "ax.tick_params(axis='both', which='major', labelsize=ticklabel_size) # Set tick label size\n",
    "plt.savefig(f'rate_dbscan_kdistance.png', dpi=300)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 1000
    },
    "id": "644iAh63gDWN",
    "outputId": "dcf7d64e-edcc-451f-e950-ca429d95918a"
   },
   "outputs": [],
   "source": [
    "# Take k-Distance around the elbow place\n",
    "elbow_index = np.argmax(np.diff(np.diff(k_distance)))  # Find the elbow point automatically\n",
    "k_distance_elbow = k_distance[elbow_index]\n",
    "print(f\"k-Distance at elbow point: {k_distance_elbow:.3f}\")\n",
    "\n",
    "# Define a range of epsilon values around the elbow point\n",
    "epsilon_range = [0.15, 0.2, 0.25]  # Use the specific values mentioned in the loop below\n",
    "print(\"Epsilon values to be used:\")\n",
    "for eps in epsilon_range:\n",
    "    print(f\"  {eps:.3f}\")\n",
    "\n",
    "for e in [0.15, 0.2, 0.25]:\n",
    "\n",
    "    # Define DBSCAN clustering model\n",
    "    mdl_dbscan = DBSCAN(eps=e, min_samples=min_pts)\n",
    "\n",
    "    # Train k-Means model\n",
    "    start_time = time.time()\n",
    "    mdl_dbscan.fit(x)\n",
    "    end_time = time.time()\n",
    "\n",
    "    print(f'Training time: {end_time - start_time:.2f} seconds')\n",
    "\n",
    "    labels_dbscan = mdl_dbscan.labels_\n",
    "\n",
    "    # Plot scatter of average rates\n",
    "    fig, ax = plt.subplots(figsize=(6,6))\n",
    "\n",
    "    # Draw noise points\n",
    "    noise_idx = (labels_dbscan == -1)\n",
    "    ax.scatter(x[noise_idx,0], x[noise_idx,1], marker=\"o\", c=\"k\", s=10**2, edgecolor=\"k\", zorder=0)\n",
    "\n",
    "    # Draw other points\n",
    "    ax.scatter(x[~noise_idx,0], x[~noise_idx,1], marker=\"o\", c=labels_dbscan[~noise_idx], s=10**2, edgecolor=\"k\")\n",
    "\n",
    "    ax.tick_params(axis='both', which='major', labelsize=ticklabel_size) # Set tick label size\n",
    "    ax.set_xlabel(f'Average rate of {genres_1}', fontsize=label_size)\n",
    "    ax.set_ylabel(f'Average rate of {genres_2}', fontsize=label_size)\n",
    "    ax.set_xlim([-0.1, 5.1])\n",
    "    ax.set_ylim([-0.1, 5.1])\n",
    "    ax.set_title(f'DBSCAN clustering, eps = {e}', fontsize=label_size)\n",
    "    ax.tick_params(axis='both', which='major', labelsize=ticklabel_size) # Set tick label size\n",
    "    plt.savefig(f'rate_dbscan_{e}.png', dpi=300) # Make figure clearer\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QNqkoRaCy80z"
   },
   "source": [
    "进一步分析最大的聚类中的用户倾向：为这些用户推荐电影"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 902
    },
    "id": "AOOxlGcUjNkd",
    "outputId": "24bd1254-b437-48cc-9cae-120ae3f3628b"
   },
   "outputs": [],
   "source": [
    "eps = 0.25\n",
    "\n",
    "# Define DBSCAN clustering model\n",
    "mdl_dbscan = DBSCAN(eps=eps, min_samples=min_pts)\n",
    "\n",
    "# Train k-Means model\n",
    "mdl_dbscan.fit(x)\n",
    "\n",
    "labels_dbscan = mdl_dbscan.labels_\n",
    "\n",
    "# Get the biggest cluster from labels_dbscan\n",
    "# Filter out noise points (labels = -1) before using np.bincount()\n",
    "biggest_cluster_label = np.argmax(np.bincount(labels_dbscan[labels_dbscan != -1]))\n",
    "print(f'Biggest cluster label: {biggest_cluster_label}')\n",
    "\n",
    "# Plot scatter of average rates\n",
    "fig, ax = plt.subplots(figsize=(6,6))\n",
    "\n",
    "# Draw noise points\n",
    "point_idx = (labels_dbscan == biggest_cluster_label)\n",
    "ax.scatter(x[point_idx,0], x[point_idx,1], marker=\"o\", c=\"k\", s=10**2, edgecolor=\"k\", zorder=0)\n",
    "\n",
    "ax.tick_params(axis='both', which='major', labelsize=ticklabel_size) # Set tick label size\n",
    "ax.set_xlabel(f'Average rate of {genres_1}', fontsize=label_size)\n",
    "ax.set_ylabel(f'Average rate of {genres_2}', fontsize=label_size)\n",
    "ax.set_xlim([-0.1, 5.1])\n",
    "ax.set_ylim([-0.1, 5.1])\n",
    "ax.set_title(f'DBSCAN clustering {biggest_cluster_label}, eps = {eps}', fontsize=label_size)\n",
    "ax.tick_params(axis='both', which='major', labelsize=ticklabel_size) # Set tick label size\n",
    "plt.savefig(f'rate_dbscan_cluster_{biggest_cluster_label}.png', dpi=300) # Make figure clearer\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "50e__pyHmuKT"
   },
   "source": [
    "用热力图的形式展示Cluster 0中用户的打分情况"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 719
    },
    "id": "W6rfuvVtmfmZ",
    "outputId": "2c7c77e2-8253-49b2-8ef7-4edadcd0e619"
   },
   "outputs": [],
   "source": [
    "# Get rate of users in cluster biggest_cluster_label\n",
    "user_ratings = ratings[ratings['userId'].isin(idx_both[labels_dbscan == biggest_cluster_label])]\n",
    "user_movie_ratings = pd.merge(user_ratings, movies, on='movieId')\n",
    "\n",
    "# Get movie list from user_movie_rating\n",
    "movie_list = user_movie_ratings['title'].unique()\n",
    "\n",
    "# Count user number grouped by movie\n",
    "user_movie_count = user_movie_ratings.groupby('title')['userId'].nunique()\n",
    "\n",
    "# Find top 20 movies names with the highest counts\n",
    "top_movies = user_movie_count.sort_values(ascending=False).head(50)\n",
    "movie_list = top_movies.index.tolist()\n",
    "\n",
    "# Filter user_movie_ratings\n",
    "user_top_movie_ratings = user_movie_ratings[user_movie_ratings['title'].isin(movie_list)]\n",
    "\n",
    "# Using heat map to show movie rate, x-axis is user_id, y-axis is movie\n",
    "user_movie_ratings_pivot = user_top_movie_ratings.pivot_table(index='userId', columns='title', values='rating', aggfunc='mean') # Change to table\n",
    "user_movie_ratings_pivot = user_movie_ratings_pivot.fillna(0) # Fill null cell with -1\n",
    "\n",
    "# Count movie number that each user has rated: rate > 0 means rated\n",
    "user_movie_count_col = user_movie_ratings_pivot.apply(lambda x: (x > 0).sum(), axis=1)\n",
    "\n",
    "# Add user_movie_count_col to user_movie_ratings_pivot as the second column\n",
    "user_movie_ratings_pivot.insert(0, 'user_movie_count', user_movie_count_col)\n",
    "\n",
    "display(user_movie_ratings_pivot)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 632
    },
    "id": "08muQng19oGY",
    "outputId": "226bc757-ce06-43ab-84ef-820568161def"
   },
   "outputs": [],
   "source": [
    "# Total user number\n",
    "print(f'Total user number: {len(user_movie_ratings_pivot)}')\n",
    "disp_user_num_left = np.floor(len(user_movie_ratings_pivot) / 2).astype(int)\n",
    "\n",
    "# Draw head image (left)\n",
    "fig, ax = plt.subplots(figsize=(100,10))\n",
    "\n",
    "im = ax.imshow(user_movie_ratings_pivot.iloc[:disp_user_num_left,1:].values.T, cmap='coolwarm') # Draw heat map\n",
    "\n",
    "# Set x-axis\n",
    "ax.set_xticks(np.arange(0, disp_user_num_left, 10))\n",
    "ax.set_xticklabels(np.arange(0, disp_user_num_left, 10))\n",
    "\n",
    "# Set y-axis\n",
    "ax.set_yticks(np.arange(len(movie_list)))\n",
    "ax.set_yticklabels(movie_list, rotation=0)\n",
    "\n",
    "# Set title and font\n",
    "ax.set_title('Heat Map of User Movie Ratings', fontsize=label_size)\n",
    "ax.tick_params(axis='both', which='major', labelsize=ticklabel_size)\n",
    "\n",
    "plt.savefig(f'rate_heatmap_left.png', dpi=300) # Make figure clearer\n",
    "plt.show()\n",
    "\n",
    "# Draw head image (right)\n",
    "fig, ax = plt.subplots(figsize=(100,10))\n",
    "\n",
    "im = ax.imshow(user_movie_ratings_pivot.iloc[disp_user_num_left:,1:].values.T, cmap='coolwarm') # Draw heat map\n",
    "\n",
    "# Set x-axis\n",
    "ax.set_xticks(np.arange(0, len(user_movie_ratings_pivot)-disp_user_num_left, 10))\n",
    "ax.set_xticklabels(np.arange(disp_user_num_left, len(user_movie_ratings_pivot), 10))\n",
    "\n",
    "# Set y-axis\n",
    "ax.set_yticks(np.arange(len(movie_list)))\n",
    "ax.set_yticklabels(movie_list, rotation=0)\n",
    "\n",
    "# Set title and font\n",
    "ax.set_title('Heat Map of User Movie Ratings', fontsize=label_size)\n",
    "ax.tick_params(axis='both', which='major', labelsize=ticklabel_size)\n",
    "\n",
    "plt.savefig(f'rate_heatmap_right.png', dpi=300) # Make figure clearer\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "按照观看电影数量进行降序排序"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 615
    },
    "id": "g-xU4VZs9pIl",
    "outputId": "3241802c-4c64-456d-d2bc-ed4cbc7f9e71"
   },
   "outputs": [],
   "source": [
    "# Sort table\n",
    "user_movie_ratings_pivot_sort = user_movie_ratings_pivot.sort_values(by='user_movie_count', ascending=False)\n",
    "user_movie_ratings_pivot_userid = user_movie_ratings_pivot_sort.index\n",
    "\n",
    "# Total user number\n",
    "disp_user_num_left = np.floor(len(user_movie_ratings_pivot_sort) / 2).astype(int)\n",
    "\n",
    "# Draw head image (left)\n",
    "fig, ax = plt.subplots(figsize=(100,10))\n",
    "\n",
    "im = ax.imshow(user_movie_ratings_pivot_sort.iloc[:disp_user_num_left,1:].values.T, cmap='coolwarm') # Draw heat map\n",
    "\n",
    "# Set x-axis\n",
    "ax.set_xticks(np.arange(0, disp_user_num_left, 10))\n",
    "ax.set_xticklabels(user_movie_ratings_pivot_userid[np.arange(0, disp_user_num_left, 10)])\n",
    "\n",
    "# Set y-axis\n",
    "ax.set_yticks(np.arange(len(movie_list)))\n",
    "ax.set_yticklabels(movie_list, rotation=0)\n",
    "\n",
    "# Set title and font\n",
    "ax.set_title('Heat Map of User Movie Ratings', fontsize=label_size)\n",
    "ax.tick_params(axis='both', which='major', labelsize=ticklabel_size)\n",
    "\n",
    "plt.savefig(f'rate_heatmap_sort_left.png', dpi=300) # Make figure clearer\n",
    "plt.show()\n",
    "\n",
    "# Draw head image (right)\n",
    "fig, ax = plt.subplots(figsize=(100,10))\n",
    "\n",
    "im = ax.imshow(user_movie_ratings_pivot_sort.iloc[disp_user_num_left:,1:].values.T, cmap='coolwarm') # Draw heat map\n",
    "\n",
    "# Set x-axis\n",
    "ax.set_xticks(np.arange(0, len(user_movie_ratings_pivot_sort)-disp_user_num_left, 10))\n",
    "ax.set_xticklabels(np.arange(disp_user_num_left, len(user_movie_ratings_pivot_sort), 10))\n",
    "\n",
    "# Set y-axis\n",
    "ax.set_yticks(np.arange(len(movie_list)))\n",
    "ax.set_yticklabels(movie_list, rotation=0)\n",
    "\n",
    "# Set title and font\n",
    "ax.set_title('Heat Map of User Movie Ratings', fontsize=label_size)\n",
    "ax.tick_params(axis='both', which='major', labelsize=ticklabel_size)\n",
    "\n",
    "plt.savefig(f'rate_heatmap_sort_right.png', dpi=300) # Make figure clearer\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "F2a1_8_4H4ez"
   },
   "source": [
    "为前用户推荐电影，即对未评分的电影预测用户可能的评分"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 525
    },
    "id": "yV00sTO3IFaM",
    "outputId": "01df78aa-d1fa-4cfb-b2bb-d977a67f9c8d"
   },
   "outputs": [],
   "source": [
    "# Select the first 100 user from user_movie_ratings_pivot_sort\n",
    "user_movie_ratings_pivot_sort_100 = user_movie_ratings_pivot_sort.iloc[:100]\n",
    "\n",
    "user_movie_ratings_pivot_userid = user_movie_ratings_pivot_sort_100.index\n",
    "\n",
    "# Draw head image (left)\n",
    "fig, ax = plt.subplots(figsize=(100,10))\n",
    "\n",
    "im = ax.imshow(user_movie_ratings_pivot_sort_100.iloc[:,1:].values.T, cmap='coolwarm') # Draw heat map\n",
    "\n",
    "# Set x-axis\n",
    "ax.set_xticks(np.arange(0, len(user_movie_ratings_pivot_sort_100), 10))\n",
    "ax.set_xticklabels(user_movie_ratings_pivot_userid[np.arange(0, len(user_movie_ratings_pivot_sort_100), 10)])\n",
    "\n",
    "# Set y-axis\n",
    "ax.set_yticks(np.arange(len(movie_list)))\n",
    "ax.set_yticklabels(movie_list, rotation=0)\n",
    "\n",
    "# Set title and font\n",
    "ax.set_title('Heat Map of User Movie Ratings', fontsize=label_size)\n",
    "ax.tick_params(axis='both', which='major', labelsize=ticklabel_size)\n",
    "\n",
    "plt.savefig(f'rate_heatmap_sort_top100.png', dpi=300) # Make figure clearer\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "计算用户相似度(Jaccard)矩阵"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "id": "q6eQcDY4I1hM"
   },
   "outputs": [],
   "source": [
    "# Compute similarity between users\n",
    "movie_list = list(user_movie_ratings_pivot_sort_100.columns[1:])\n",
    "user_list = list(user_movie_ratings_pivot_sort_100.index)\n",
    "sim_matrix = np.zeros((len(user_list), len(user_list)))\n",
    "\n",
    "for i in range(len(user_list)):\n",
    "    for j in range(len(user_list)):\n",
    "        if i == j: # Avoid the influence of itself\n",
    "            continue\n",
    "        rate_i = user_movie_ratings_pivot_sort_100.iloc[i,1:]\n",
    "        rate_j = user_movie_ratings_pivot_sort_100.iloc[j,1:]\n",
    "        intersection = np.sum((rate_i > 0) & (rate_j > 0))\n",
    "        union = np.sum((rate_i > 0) | (rate_j > 0))\n",
    "        sim_matrix[i,j] = intersection / union if union > 0 else 0\n",
    "\n",
    "# Get user_list_unrate to store users who have movies not been rated\n",
    "user_list_unrate = {}\n",
    "for i in range(len(user_movie_ratings_pivot_sort_100)):\n",
    "    movie_unrate = np.where(user_movie_ratings_pivot_sort_100.iloc[i,1:] == 0)[0] # Start from the first movie while jumping over the first column\n",
    "    if len(movie_unrate) > 0:\n",
    "        user_id_unrate = user_movie_ratings_pivot_sort_100.index[i]\n",
    "        user_list_unrate[user_id_unrate] = np.where(user_movie_ratings_pivot_sort_100.iloc[i,1:] == 0)[0] + 1 # Add 1 to real pos\n",
    "\n",
    "# Print rate of unrated movies\n",
    "# for unrate_user_id in user_list_unrate.keys():\n",
    "#     for unrate_movie_id in user_list_unrate[user_id]:\n",
    "#         print(user_movie_ratings_pivot_sort_100.iloc[user_list.index(user_id), unrate_movie_id])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "gcC_CJ6Ym--l"
   },
   "source": [
    "对单个例子进行评分"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "id": "7rmiCZETdSyE",
    "outputId": "fb68b6da-3924-4e89-d50b-bdf0cdc942e5"
   },
   "outputs": [],
   "source": [
    "def rate_movie(unrate_user_id, unrate_movie_id, user_movie_ratings_pivot_sort_100, sim_matrix):\n",
    "    co_user_num = 10 # Colabarative users - k value\n",
    "\n",
    "    # User idx who rate unrate_movie_id[i] higher than 0\n",
    "    unrate_movie_rates = user_movie_ratings_pivot_sort_100.iloc[:,unrate_movie_id]\n",
    "    co_user_idx = np.where(unrate_movie_rates > 0)[0]\n",
    "\n",
    "    # Get co-users' rate on the unrated movie\n",
    "    co_user_rate = unrate_movie_rates.iloc[co_user_idx]\n",
    "\n",
    "    # Get co-users' similarity to the unrate_user_id\n",
    "    co_user_sim = sim_matrix[user_list.index(unrate_user_id), co_user_idx]\n",
    "\n",
    "    # Create co_user_rate_sim by combining co_user_rate with co_user_sim\n",
    "    co_user_rate_sim = np.array([co_user_rate, co_user_sim]).T\n",
    "\n",
    "    # Sort co_user_rate_sim\n",
    "    co_user_rate_sim_sort = co_user_rate_sim[co_user_rate_sim[:,1].argsort()[::-1]]\n",
    "\n",
    "    # Predict movie rate\n",
    "    top_rate = co_user_rate_sim_sort[:co_user_num,0]\n",
    "    top_sim = co_user_rate_sim_sort[:co_user_num,1]\n",
    "    predict_rate = np.dot(top_rate, top_sim) / np.sum(top_sim)\n",
    "\n",
    "    return predict_rate\n",
    "# Single list\n",
    "unrate_user_id = list(user_list_unrate.keys())[0]\n",
    "unrate_movie_id = user_list_unrate[unrate_user_id][0]\n",
    "unrate_movie_name = movie_list[unrate_movie_id]\n",
    "print(f'User {unrate_user_id} has {unrate_movie_name} ({unrate_movie_id}) unrated movies')\n",
    "\n",
    "predict_rate = rate_movie(unrate_user_id, unrate_movie_id, user_movie_ratings_pivot_sort_100, sim_matrix)\n",
    "print(f'Predicted rate of {unrate_movie_name} ({unrate_movie_id}) for user {unrate_user_id} is {predict_rate}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "EdzmRwPmnGhk"
   },
   "source": [
    "对所有未评分用户进行电影喜好预测"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 525
    },
    "id": "hPLx40LnnOei",
    "outputId": "109e0f6c-0d90-42d4-ca68-13b4fb5514d6"
   },
   "outputs": [],
   "source": [
    "user_movie_ratings_pivot_sort_100_fill = user_movie_ratings_pivot_sort_100.copy()\n",
    "\n",
    "# Predict movie rate of all unrated movies\n",
    "for unrate_user_id in user_list_unrate.keys():\n",
    "    for unrate_movie_id in user_list_unrate[unrate_user_id]:\n",
    "        rate_pred = rate_movie(unrate_user_id, unrate_movie_id, user_movie_ratings_pivot_sort_100, sim_matrix)\n",
    "\n",
    "        # Fill user_movie_ratings_pivot_sort_100\n",
    "        user_movie_ratings_pivot_sort_100_fill.iloc[user_list.index(unrate_user_id), unrate_movie_id] = rate_pred\n",
    "\n",
    "# Draw head image (left)\n",
    "fig, ax = plt.subplots(figsize=(100,10))\n",
    "\n",
    "im = ax.imshow(user_movie_ratings_pivot_sort_100_fill.iloc[:,1:].values.T, cmap='coolwarm') # Draw heat map\n",
    "\n",
    "# Set x-axis\n",
    "ax.set_xticks(np.arange(0, len(user_movie_ratings_pivot_sort_100_fill), 10))\n",
    "ax.set_xticklabels(user_movie_ratings_pivot_userid[np.arange(0, len(user_movie_ratings_pivot_sort_100_fill), 10)])\n",
    "\n",
    "# Set y-axis\n",
    "ax.set_yticks(np.arange(len(movie_list)))\n",
    "ax.set_yticklabels(movie_list, rotation=0)\n",
    "\n",
    "# Set title and font\n",
    "ax.set_title('Heat Map of User Movie Ratings', fontsize=label_size)\n",
    "ax.tick_params(axis='both', which='major', labelsize=ticklabel_size)\n",
    "\n",
    "plt.savefig(f'rate_heatmap_sort_top100_fill.png', dpi=300) # Make figure clearer\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 实验七"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " 实现k-Means聚类算的核心功能，并使用虚拟数据进行验证\n",
    " 使用肘部法确定k值\n",
    "使用OpenCV库读取图像数据，并转换为灰度图像\n",
    "使用k-Means生成蒙版，区分灰度图像的前景与背景\n",
    "使用k-Means对像素颜色进行聚类，观察不同k值条件下的结果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import cv2\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.cluster import KMeans\n",
    "from sklearn.metrics import pairwise_distances_argmin_min\n",
    "\n",
    "# 1. 实现k-Means聚类算法的核心功能\n",
    "\n",
    "def k_means(X, k, max_iter=100, tol=1e-4):\n",
    "    \"\"\"\n",
    "    实现k-Means算法。\n",
    "    :param X: 输入数据，二维数组，形状为(n_samples, n_features)\n",
    "    :param k: 聚类中心的数量\n",
    "    :param max_iter: 最大迭代次数\n",
    "    :param tol: 收敛容忍度\n",
    "    :return: labels (聚类标签), centers (聚类中心)\n",
    "    \"\"\"\n",
    "    # 随机初始化聚类中心\n",
    "    n_samples, n_features = X.shape\n",
    "    centers = X[np.random.choice(n_samples, k, replace=False)]\n",
    "    \n",
    "    for i in range(max_iter):\n",
    "        # 计算每个样本点到聚类中心的距离，并分配标签\n",
    "        labels = pairwise_distances_argmin_min(X, centers)[0]\n",
    "        \n",
    "        # 更新聚类中心\n",
    "        new_centers = np.array([X[labels == j].mean(axis=0) for j in range(k)])\n",
    "        \n",
    "        # 判断是否收敛\n",
    "        if np.linalg.norm(new_centers - centers) < tol:\n",
    "            break\n",
    "        centers = new_centers\n",
    "    \n",
    "    return labels, centers\n",
    "\n",
    "# 验证k-Means聚类算法（使用虚拟数据）\n",
    "from sklearn.datasets import make_blobs\n",
    "\n",
    "X, y = make_blobs(n_samples=300, centers=4, random_state=42)\n",
    "\n",
    "# 使用k-Means聚类\n",
    "k = 4\n",
    "labels, centers = k_means(X, k)\n",
    "\n",
    "# 绘制聚类结果\n",
    "plt.scatter(X[:, 0], X[:, 1], c=labels, cmap='viridis')\n",
    "plt.scatter(centers[:, 0], centers[:, 1], c='red', marker='x', s=200, label='Centers')\n",
    "plt.title(\"k-Means Clustering (Virtual Data)\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "# 2. 使用k-Means进行图像的背景分割\n",
    "\n",
    "def k_means_image_segmentation(image_path, k, max_iter=100):\n",
    "    \"\"\"\n",
    "    使用k-Means进行图像的背景分割。\n",
    "    :param image_path: 图像路径\n",
    "    :param k: 聚类数量\n",
    "    :param max_iter: 最大迭代次数\n",
    "    :return: 分割后的图像\n",
    "    \"\"\"\n",
    "    # 读取图像并转换为RGB格式\n",
    "    image = cv2.imread(image_path)\n",
    "    image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)\n",
    "    \n",
    "    # 将图像转换为二维数据，每一行代表一个像素的RGB值\n",
    "    pixels = image.reshape(-1, 3)\n",
    "    \n",
    "    # 执行k-Means聚类\n",
    "    labels, centers = k_means(pixels, k, max_iter)\n",
    "    \n",
    "    # 根据聚类结果重新构建图像\n",
    "    segmented_image = centers[labels].reshape(image.shape)\n",
    "    \n",
    "    return segmented_image.astype(np.uint8), labels.reshape(image.shape[:-1])\n",
    "\n",
    "# 读取并分割图像\n",
    "image_path = 'path_to_image.jpg'  # 请替换为图像路径\n",
    "k = 2  # 选择前景和背景的聚类数量\n",
    "segmented_image, mask = k_means_image_segmentation(image_path, k)\n",
    "\n",
    "# 显示分割结果\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.imshow(segmented_image)\n",
    "plt.title(f\"Segmented Image (k={k})\")\n",
    "plt.axis('off')\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.imshow(mask, cmap='gray')\n",
    "plt.title(\"Mask (Foreground vs Background)\")\n",
    "plt.axis('off')\n",
    "\n",
    "plt.show()\n",
    "\n",
    "# 3. 使用k-Means进行图像压缩\n",
    "\n",
    "def k_means_image_compression(image_path, k, max_iter=100):\n",
    "    \"\"\"\n",
    "    使用k-Means进行图像压缩。\n",
    "    :param image_path: 图像路径\n",
    "    :param k: 聚类数量（压缩等级）\n",
    "    :param max_iter: 最大迭代次数\n",
    "    :return: 压缩后的图像\n",
    "    \"\"\"\n",
    "    # 读取图像并转换为RGB格式\n",
    "    image = cv2.imread(image_path)\n",
    "    image = cv2.cvtColor(image, cv2.COLOR_BGR2RGB)\n",
    "    \n",
    "    # 将图像转换为二维数据，每一行代表一个像素的RGB值\n",
    "    pixels = image.reshape(-1, 3)\n",
    "    \n",
    "    # 执行k-Means聚类\n",
    "    labels, centers = k_means(pixels, k, max_iter)\n",
    "    \n",
    "    # 根据聚类结果重新构建图像\n",
    "    compressed_image = centers[labels].reshape(image.shape)\n",
    "    \n",
    "    return compressed_image.astype(np.uint8)\n",
    "\n",
    "# 读取并压缩图像\n",
    "compressed_image = k_means_image_compression(image_path, k=8)\n",
    "\n",
    "# 显示压缩结果\n",
    "plt.figure(figsize=(10, 5))\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.imshow(cv2.imread(image_path))\n",
    "plt.title(\"Original Image\")\n",
    "plt.axis('off')\n",
    "\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.imshow(compressed_image)\n",
    "plt.title(f\"Compressed Image (k=8)\")\n",
    "plt.axis('off')\n",
    "\n",
    "plt.show()\n"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "authorship_tag": "ABX9TyO+9EvQyLwzW6sRWnzqnjPu",
   "collapsed_sections": [
    "c-Xj4WtjPdwD"
   ],
   "include_colab_link": true,
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
